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A short gas discharge layer sandwiched with a semiconductor layer between planar electrodes 
shows a variety of spatio-temporal patterns. The paper focusses on the spatially homogeneous 
spontaneous oscillations while a DC voltage is applied; the results on these homogeneous oscillations 
apply equally to a planar discharge in series with any resistor with capacitance. We define the 
minimal model, identify its independent dimensionless parameters and then present results of the 
full time-dependent numerical solutions of the model as well as of a linear stability analysis of the 
stationary state. Full numerical solutions and the results of the stability analysis agree very well. 
The stability analysis is then used for calculating bifurcation diagrams. We find semi-quantitative 
agreement with experiment for the diagram of bifurcations from stationary to oscillating solutions 
as well as for amplitude and frequency of the developing limit cycle oscillations. 



I. INTRODUCTION 

Gas discharges on the transition from Townsend to 
glow regime exhibit a wealth of spatio-temporal struc- 
tures. Besides striations, i.e^longitudinal waves in a long 
discharge column 0,0,|3)S@1j short discharges with wide 
lateral aspect ratio can also exhibit rich spatio-temporal 
structures in the transversal direction as reported by a 
number of authors 0, 0, H, I3- This is even the case 
when the externally applied voltage is stationary and the 
gas is pure, as long as the system is sandwiched between 
planar electrodes and at least one Ohmic layer. An in- 
teresting seq uence of experiments has been performed in 
Miinster where the bifurcations between different 

spatio-temporal states in parameter space were investi- 
gated very systematically. 

As in our previous paper |l2l |. we focus in the present 
one on the purely temporal oscillations that occur in a 
spatially homogeneous mode. This focus has two reasons: 
first, understanding the temporal structures is a first 
systematic step towards understanding the full spatio- 
temporal structures; second, there are numerous obser- 
vations of temporal oscillations in comparable parameter 
regimes [H El 111 E [H El El ^ . For the oscilla- 
tions, the setup need not contain an Ohmic layer as in 
[Tol a resistor with capacitance in the circuit will 
have the same effect on the gas discharge. 

In the previous paper |l2l |. we concentrated on the 
question whether a simple two-component reaction- 
diffusion model for current and voltage in the gas dis- 
charge layer would be sufficient to describe the oscilla- 
tions. Such a model is suggested through similarities 
with patterns formed in a number of physical, chemical 
or biological systems like the Belousov-Zhabotinski reac- 
tion, Rayleigh-Bcnard convection, patterns in bacterial 
colonies, in Dictyostelium or in nerval tissue etc. How- 
ever, the actual results of a realistic gas discharge model 
are in conflict with a simple two-component reaction dif- 



fusion approximation that neglects the height and subse- 
quent memory of the system. This can be seen, in par- 
ticular, from the occurence of a period dou bling cascade 
as well as from analytical model reductions Tz| . Similar 
period doublin g ca scades are observed expcrimcntaly in 

lUilliKililH. 

In the present paper, we continue the analysis of the 
full gas discharge model, coupled to a high-Ohmic layer 
and driven by a stationary voltage. The focus is now on 
quantitative comparison with experiment, on a stability 
analysis and on the derivation of bifurcation diagram. 
The specific experiment to be analyzed was performed in 
nitrogen at 40 mbar within a gap of 0.5 or 1 mm wide 
while the semiconductor was a layer of 1.5 mm photosen- 
sitively doped GaAs. To the whole structure, voltages in 
the ran ge o f 500 to 800 V were applied. As in our previous 
papers 123, j restrict the analysis to the direc- 
tion normal to the layers, hence assuming homogeneity 
in the transversal directions. The experimental system 
actually shows a transition from a homogeneous station- 
ary to a homogeneous oscillating state, and the theory 
presented here reproduces essential features of these ex- 
periments. At the same time, the investigation serves 
as a gauge point for a later analysis of spatio-temporal 
patterns. 

In detail, we define the model as a set of partial differ- 
ential equations and perform a dimensional analysis in 
Section II. In Section III, first the physical parameters 
and the numerical details of solving the PDE's in time 
are given. Then qualitative and quantitative results of 
numerical solutions and experiments are discussed. In 
particular, the hysteresis between stationary and oscil- 
lating solutions is demonstrated numerically, amplitude 
and frequency of the limit cycle oscillations as a function 
of applied voltage and conductivity of the semiconductor 
are compared with experimental results, and the physi- 
cal mechanism of the oscillation is discussed. In Section 
IV, it is explained how the stability analysis about a sta- 



2 



tionary solution of the complete system is performed. In 
Section V, the results of the stability analysis are pre- 
sented. First a convincing agreement between numerical 
solutions of the full PDE's and the stability analysis re- 
sults is found. Then the stability analysis is used to cal- 
culate bifurcation diagrams for the transition from sta- 
tionary to oscillating states that are then compared with 
experiment. The paper concludes with Section VI. 

II. THE MODEL 

The experiment consists of two layers, a gas discharge 
and a semiconductor, sandwiched between two planar 
electrodes to which a DC voltage is applied. In this sec- 
tion, the equations are defined and a dimensional anal- 
ysis is performed to identify the independent parameter 
combinations of the problem. This also serves to identify 
physical processes and time scales. 

A. Gas discharge layer 

In the gas discharge, two ionization mechanisms coop- 
erate to maintain conductivity: the so-called a process of 
impact ionization in the bulk of the discharge, and the 7 
process of secondary emission at the cathode. The clas- 
sical "fluid" approximation consists of continuity equa- 
tions for electron density Uf, and positive ion density n+ , 
coupled to the Poisson equation for the electric field E: 



dt Tie -f drJe — SOUrCC , (1) 

dt n+ + drJ+ = source , (2) 
drE = — (n+ - ne) . (3) 



The spatial coordinate r is normal to the layers, and in 
the present paper, it is assumed that there are no vari- 
ations in the transversal directions. The gas is assumed 
to be non-attaching, i.e., no negative ions are formed. 
Also photo-ionization, Ohmic heating, nonlocal interac- 
tions and diffusion are neglected in this simplest approx- 
imation. The particle current densities Je and J+ are 
approximated by a drift motion that is linear in the field 

Je = -He fie E , J+ ^ n+ /i+ E . (4) 

The source term on the right hand side of Eqs. and 
(O is approximated by impact ionization in the classical 
Townsend form 

source = |neMe£^| e"-^"/'^' . (5) 

The one-dimensional approximation of Eqs. (QJ, Q 

and Q makes the total electric current J{t) homoge- 
neous 

endtE{r,t) +eJe{r,t) +eJ+{r,t) ^ J{t) , drJ{t) = 0. 

(6) 



This identity can be used to substitute Je or J+ by J(t). 
In the present analysis, we will keep ne{r,t) and E{r,t) 
as independent fields and express n+ and J_|_ by these 
fields and the total current J{t). 

The model is completed by boundary conditions on 
the electrode. At the anode which is located at r = 0, 
electrons are absorbed and ions are absent: 

J+(0,i) = <^ n+(0,t) = 0. (7) 

At the cathode which is located aX r = d, impacting ions 
can liberate electrons by secondary emission with rate 7: 

\Je{d,t)\ =7 I J+ (d, t) I ^e^e {d, t) = ^n+n+ (d, t) . 

(8) 

Note that consistenly with [T^ l2^ , but in contrast with 
most other literature, the anode is on the left hand side 
at r = 0. This has the advantage that the electric field 
is positive, and sign mistakes when evaluating E or \E\ 
cannot occur. 

Substantial densities of charged particles change the 
electric field according to ©, and the electric field deter- 
mines drift and ionization rates of the particles according 
to Eqs. CJ, ©I Q and (0). Therefore the process is non- 
linear as soon as space charges become relevant. It causes 
the well-known transition from the linear Townsend dis- 
charge to the nonlinear glow discharge. 



B. Semiconductor layer and complete circuit 

The semiconductor layer of thickness dg is assumed to 
have a homogeneous and field independent conductivity 
as and dielectricity constant e^: 

Jsit) ^ a^Esit) , q = eseodrE. (9) 

As there are no space charges in the bulk of the semicon- 
ductor, the electric field is homogeneous, and voltage and 
field are related through Us{t) = Es{t)ds. The equation 
of charge conservation dtq + drJs = in one dimension 
leads again to the homogeneity of the total current den- 
sity J{t) 

ese^dtEs{t) + Js{t) = J{t), (10) 

that is the same as in the gas discharge 0. Hence in 
macroscopic parameters, the semiconductor solves 

CsdtUs{t) + JS)^J{t) , C/.(t) = i?sJs(t), (11) 

= ^ , Rs^^. (12) 
ds (Js 

where Cg is the capacitance per area. 

According to ((TT)l . perturbations of Us{t) or Js[t) de- 
cay on the Maxwell time scale 

= CsRs = — . (13) 
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This time scale is independent of the thickness of the 
semiconductor layer although it represents the time that 
the charge needs to cross it. The time scale of the exper- 
imentally observed oscillations is of the order of Tj, and 
therefore also approximately proportional to 1 j Og as will 
be discussed in Section III.D. 

Actually, for the present investigation of one- 
dimensional oscillations, the specific structure of a planar 
semiconductor layer is not required, but any serial com- 
ponent of the electric circuit with capacitance Cs and 
resistance i?s will support the same equation 

The total stationary voltage Ut over the complete sys- 
tem is 

Ut = U(t)^Us{t) , U{t)^ f E{r,t)dr , dtUt^O. 



(14) 

According to and 1)14(1 . the dynamics of the voltage 
U{t) on the gas discharge obeys the equation 



TsdtU ^Ut- U{t) - R,J{t). 



(15) 



C. Dimensional analysis and system definition 

The dimensional analysis is performed as previously in 
[T^ ISS 22\ ■ We introduce the dimensionless coordi- 
nates and fields 

z = ^ , , cr{z,T)=^ , (16) 

Aq to no 



J 



Eo 



EqXo 



enoXo/to ' 



measuring quantities in terms of the intrinsic parameters 
of the system 



^0 — — , ^0 — — 



no = . (17) 



After eliminating the ion dynamics by the total current 
j (r) , the equation of motion of the gas discharge becomes 

dr<T = dzje + jea{£) , je = oE, (18) 

drS - 3{t) - (1 + p)3e - fi£d,£, (19) 



and the boundary conditions iQ and © read 



dr£{0,T) - j(T)-ie(0,r), 



(20) 



dr£{L,T) = j{t) 



1 + 7 
7 



Je{L,T). (21) 



The intrinsic dimensionless parameters of the gas dis- 
charge are the mobility ratio /i of electrons and ions and 
the length ratio L of system size and inverse cross section 
of impact ionization 



Me Ao 



(22) 



The discharge is coupled to the semiconductor and the 
DC voltage source Ut through ((TT|l as 

TsdMir) ^ Ut- U{t) - ns]{T), (23) 

with the dimensionless parameters 

Ts _ Rs 



to ' '^'^ Eoto/{cno)' 



(24) 



The voltage W(t) = £{ztt) dz is related to the electric 
field £ and potential in differential form as 

£{z,t) = -d,^{z,T) , U{t)=4>{Q,t)-4>{L,t), (25) 

where gauge freedom allows one to choose 

0(O,r)=O. (26) 

Hence the dynamics of the complete system is de- 
scribed by Eqs. dll-lEU, ^ &nd The sys- 
tem is characterized completely by the independent di- 
mensionless parameters /i, L and 7 for the gas discharge 
layer, Tj and TZg for the semiconductor layer and the total 
applied DC voltage Ut- 



III. NUMERICAL SOLUTIONS OF THE 
DYNAMICS 

In this section, this dynamical model is solved numeri- 
cally and the results are compared with experiments. We 
discuss physical parameters under A and numerical de- 
tails under B. In C, qualitative features of experimental 
and numerical system are compared like the bistability 
between stationary and oscillating state. In D, a quanti- 
tative comparison between theory and experiment is per- 
formed, and the dependence of amplitude and frequency 
of the oscillation as a function of Ut and l/TZg is deter- 
mined numerically. Finally, in E, we discuss the mech- 
anism of the oscillations and identify the surface charge 
effects that are inherent in our model. 



A. Physical parameters 

In the experiment |0, nitrogen at a pressure of 40 
mbar was used in gaps with widths of 0.5 or 1 mm. The 
article ^3 contains mainly data for the 0.5 mm gap, 
while the Ph.D. thesis [llj also contains more data for 
1 mm. The gas discharge was coupled to a semiconduc- 
tor layer of GaAs with a width of = 1.5 mm and a 
dielectricity constant Cs — 13.1. Through photosensi- 
tive doping, the conductivity of the semiconductor layer 
could be increased by about an order of magnitude; the 
dark conductivity was ag — 3.2 • 10~*(r2cm)~^. For the 
discharge gap of 0.5 mm width, voltages in the range of 
500 to 600 V were used; for the gap of 1 mm width, the 
applied voltages were in the range of 580 to 740 V. 



4 



Of course, the predictive power of the theory depends 
on the model approximations as weh as on the chosen 
parameters. Our simple classical model will not give fully 
quantitative agreement. On the other hand, its simple 
structure and few parameters give a chance of physical 
understanding and control. 

For the gas discharge, we used the ion mobility ^.j^ = 
23.33 cm^/Vs and electron mobility = 6666.6 cm^/Vs. 
For ao ^ Ap = [27.78/im]-i and Eq = Bp = 
10.26 kV/cm, the value from [s^l was used. The gap 
widths of d = 0.5 and 1 mm then correspond to dimen- 
sionless gap widths L = 18 and 36. For 7, we used the 
value 0.08 determined from experimental Paschen curves 
in 0. It should be noted that our classical model pre- 
dicts that the Paschen curves (i.e., the breakdown volt- 
age U of the gas discharge as a function of pressure times 
gap width pd) for different system sizes should be indis- 
tinguishable. In practice, they do not precisely fall on 
top of each other. 
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FIG. 1: Current-voltage characteristics for 7 = 0.08 (solid 
lines) and 7 = 0.1 (dashed line) for the dimensionless gap 
widths L as indicated in the figure. 



voltage has only one minimum as function of current j 
and this occurs for a value j ^ 0. This subcritical be- 
havior corresponds to the classical textbook case where 
the characteristics bends down from the Townsend break- 
down voltage towards a voltage minimum in th e gl ow 
discharge regime — as we have discussed in [23, IHj in 
detail, this requires a sufficiently large system size. For 
7 = 0.08, the characteristics becomes subcritical for sys- 
tem size L > Lcrit = In [(1 -I- 7)/7] — 19.2 while the 
transition to supercritical behavior is determined numer- 
ically I23 to the value of L = 17.2. 

Data on the coefficient 7 of secondary electron emis- 
sion are relatively scarce, so it is quite common 31] to 
use it as an adjustable parameter as we do. The tabu- 
lated data for ao = Ap and Eq = Bp from 13(1 to gether 
with the Paschen curve for d = 0.5 mm from ll'l would 
suggest 7 — 0.03, but that would mean that the char- 
acteristics would be supercritical up to i = 24.9, then 
it would develop some regime with negative differential 
conductivity, and it would become subcritical only for 
L > Lc„t{l = 0.03) = 26.1. 

We conclude that the gap with width 0.5 mm (cor- 
responding to L = 18) is so sensitive to the not very 
well known parameter 7 that an analysis of the exper- 
imental data would be rather uncertain. Furthermore, 
the approximation of purely local interactions becomes 
worse in shorter gaps. Finally, the electric fields in short 
discharges are higher and vary more; therefore the as- 
sumption that 7 does not depend on E becomes more 
restrictive. For this reason, we chose to analyze the sys- 
tem with gap width 1 mm (L = 36). 

We recall that the following intrinsic scales 

Xa « 27.78 ^m , to- 40.6 • lO^^^ 
no ~ 2.04 • lO^Vcm^ , Ea w 10.26 kV/cm (27) 

enter the dimensional analysis H16|l . Therefore the di- 
mensionless parameters for a system with gap width of 
d = 1 mm and applied voltages in the range from 500 to 
740 V are in our simulations: 



It is interesting to note how sensitive the theoretical 
results are to small changes of the secondary emission co- 
efficient 7, in particular, for the short gap with L — 18. 
This is illustrated in Fig. 1. The upper three solid lines 
show the shape of the current voltage characteristics for 
7 = 0.08 and gap widths of L = 17, 17.5 and 18. As dis- 
cussed in more detail in [23, 13 1 ^^e characteristics can 
be supercritical {L — 17, positive differential conductiv- 
ity for all values of the current j), mixed II (L = 17.5, 
Townsend breakdown voltage lower than the local volt- 
age minimum for j ^ 0) or mixed I (L = 18, Townsend 
breakdown voltage higher than the local voltage mini- 
mum for j ^ 0). The dashed line shows the characteris- 
tics for L =18 and 7 = 0.1. U then overall is considerably 
lower and the characteristics is fully subcritical, i.e., the 



Ai = 0.0035 , L = 36 , 7 = 0.08, 

T, = 0.243 TZs , 3 • 10^ < 7^s < 3 • 10^ 

17.5 <Ut< 26. (28) 

Here, the dimensionless capacitance of the semiconductor 
layer is Cg — 0.243, and its dimensionless characteristic 
time scale is Ts = CsTZs- The value TZg = 3 ■ 10^ for the 
semiconductor resistance corresponds to the dark con- 
ductivity of as = 3.2 • 10-s/(f]cm), and 7^s = 3 • 10^ 
corresponds to the fully photo-activated conductivity 
(Ts = 3.2 • 10~^/(ricm). The dimensionless voltage range 
of 17.5 <Ut<26 corresponds to the dimensional range 
of 500 Y <Ut< 740 V. 
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B. Numerical solution strategy 

Equations lfTS|) - (PH|l were solved numerically with an 
implicit temporal discretization, which makes the calcu- 
lation numerically stable for arbitrary time and space 
steps. After discretization, the dynamical equations H18f) 
and H19|l have the form 



m+l 
1 



m+l 



At 



At 



Az 



-1 c'Tl+l 



+ 



(29) 



where i parametrizes the spatial and m the temporal grid. 

For known ct™ and at time step m, the boundary 
condition on the left (|20|1 determines 



^1 ^ "--l 



A 



(30) 



then the other fields 5™^^ are calculated successively 
from the left to right {i = 2, 3, .., N) by the equation 



?m+l 



Sr f 1 + '^Sr^' - (1 + Ara™) + Arj' 



1 -X- ifAr Cit: 



(31) 

For cr™^^, the boundary condition on the right (|21() de- 
termines 
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"AT 



(32) 



The remaining cr™^^ can now be calculated successively 
from the right to left {i = N -1,N -2, .., 1) as 



m+l 



(33) 



1 + ^£r+^ - Arf ™+ia(fr+') ■ 

The total current in these equations is determined by 



1 



J 



TZs + TsL 



+ (l + /i)Az^(£a)r 

i=l / . 



.(34) 



This identity can be derived from (|23|l where drU is iden- 
tified with dz dr£ through H25|) . and then for dr£, the 
identity p9|) is used. 

The results presented in Figures 2 to 9 are derived on 
a grid with Az = 36/600 and At — 180/600 which gives 
a sufficient numerical accuracy. 



C. Qualitative features of experimental and 
numerical oscillations: hysteresis amd limit cycles 

The experiments [lOj show approximately periodic os- 
cillations. They are quite anharmonic with long phases 



of low current interrupted by a short current pulse. De- 
pending on applied voltage Ut and resistance of the semi- 
conductor layer TZs , either the homogeneous stationary or 
the homogeneous oscillating state are dynamically sta- 
ble. Inbetween, there is a regime of bistability where it 
depends hysteretically on the previous state whether the 
system is stationary or oscillating. 

The same qualitative behavior can be observed in our 
numerical solutions. First, the upper panel in Fig. 2 
shows the current j(r) as a function of time for the sys- 
tem with the parameters from (|28l) and TZg = 4:- 10^ and 
Ut — 19.5 (which corresponds to as = 2.4 • 10^^/(ricm) 
and Ut = 555 V). After some transient, the current 
relaxes to periodic unharmonic oscillations. The lower 
panel in Fig. 2 shows the voltage U{t) over the gas dis- 
charge; the voltage on the semiconductor is correspond- 
ingly Ut — U{t). In dimensional units, the peak current 
of the oscillations is about 9 mA/cm^ and the frequency 
is about 120 kHz. 



,x 10 




X 10 



FIG. 2; j(r) and W(r) for the parameters from I28I I. TZs 
4 • 10^ and Ut = 19.5. 



The same numerical data for current j and voltage 
U are shown as a phase space plot in Fig. 3. The fig- 
ure shows more precisely the approach to a limit cycle. 
Fig. 3 contains two additional lines, namely the current 
voltage characteristics of the gas discharge U = U{j) and 
the load line U = Ut — TZsj- Their intersection marks the 
stationary solution of the system. In the present case, it 
is located in the low current regime close to the Townsend 
limit, while the peak current explores the regime of sub- 
normal glow. 

The system of Figs. 2 and 3 is actually in the bistable 
regime. For different initial conditions that are a suffi- 
ciently small perturbation of the stationary state, the 
same system relaxes to the stationary point. This is 
shown as phase space plot in Fig. 4. 
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FIG. 3: Phase space plot of the data from Fig. 2. After some 
transient time, a stable limit cycle is reached. Also drawn 
are the current-voltage-characteristics U = U{j) of the gas 
discharge and the load line U — Ut — TZsj- Their intersection 
denotes the stationary solution. 



of a system where the stationary solution is dynamically 
unstable, and the system runs away from this initial state 
and eventually reaches a limit cycle oscillation. This be- 
havior is shown in Fig. 5 as j(r) and U{t), while Fig. 6 
shows the corresponding phase space plot. 
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FIG. 5; j(r) and W(r) for the parameters from (I28I I. TZs = 
4 • 10^ and Ut ~ 24. The stationary state now is linearly 
unstable and develops into a limit cycle. 




FIG. 4: System with exactly the same parameters as in Figs. 2 
and 3, but for different initial conditions. The system now 
spirals inwards towards the stationary point. 



If the applied voltage Ut becomes large enough, the 
stationary state becomes unstable for any initial condi- 
tion. The search for appropriate parameters was guided 
by the stability analysis described in sections IV and V 
of this paper. We find that Ut = 24 {Ut = 684 V) with all 
other parameters unchanged can be used as an example 



FIG. 6: Phase space plot of the data from Fig. 5 with current- 
voltage-characteristics and load line. 
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D. Quantitative comparison: amplitude and 
frequency of oscillations 

The qualitative agreement of numerical solutions and 
experiment now encourages a more quantitative compar- 
ison. The thesis contains diagrams on how frequency 
and maximal current amplitude depend on the semicon- 
ductor conductivity for a gas gap of 1 mm. It also con- 
tains the remark that frequency and amplitude for fixed 
conductivity depend in about the same way on the ap- 
plied voltage as in the 0.5 mm gap of Ref. |10l |. 

The same diagrams can also be derived from the nu- 
merically obtained limit cycle oscillations, they are pre- 
sented in Fig. 7. The figure shows the current amplitude 
A and frequency / as a function of semiconductor con- 
ductance 1/TZs for fixed voltage Ut or as a function of Ut 
for fixed l/TZg- 



X 10 



X 10" 





FIG. 7: Amplitude A and frequency / of the current oscilla- 
tions as a function of applied total voltage Ut (for fixed resis- 
tance Ra = 4 ■ 10^) and as a function of conductivity l/TZa 
(for fixed voltage Ut ~ 21). 



We now compare the results. The upper left panel 
shows that the maximal current amplitude A as a func- 
tion of applied voltage Ut is increasing with decreasing 
slope. This agrees with the statements written in [lH . 
The upper right panel shows that the frequency / is an 
almost linearly increasing function of applied voltage Z^, 
this is actually in contradiction with the statement in 
that the function would decrease. 

The lower two panels allow a more quantitative com- 
parison since corresponding experimental diagrams can 
be found in pdj . The experiments explore the range of 
0.6 • 10"^/(f7cm) < as < 2.8 • 10"^/(f7cm) which corre- 
sponds to 0.62- 10"^ < l/TZs < 2.9-10"^. The experimen- 
tal diagrams for Ut = 605 V and 616 V in show, that 
the amplitude A is very sensitive to this change while the 



frequency / is rather robust. The numerical results are 
derived for Z/^t = 21 which corresponds to Ut — 600 V. 

In detail, the experimental curve for the current am- 
plitude for 605 V shows first an increase from 0.2 to 0.8 
mA with a subsequent sudden drop to essentially from 
which the current suddenly jumps to values from 1.0 to 
1.5 mA. For 616 V, in contrast, an almost continuous 
increase from 0.2 to 2.7 mA is observed for the same 
resistance range. Not too suprisinngly, our numerical re- 
sults reproduce neither of these widely differing results 
at quite similar voltage. Rather, we observe an almost 
constant value in the range of 5.5 • 10""* to 6.0 • 10^"* in 
the lower left pannel. 

On the other hand, for the variation of the frequency 
/ with conductivity, experiments JXl'l both for 605 V and 
for 616 V observe an about linear increase from 115 kHz 
or 125 kHz to 220 kHz (4.6 • 10"^ < / < 8.8 • lO"*^ in 
our dimensionless units) in the range of 0.62 • 10~^ < 
l/TZs < 2.9 • 10^^. Our numerical results in this range 
of l/TZs show the same linear increase, from 1.5 • 10~® 
to 6.5 • 10~^. We believe that this agreement is quite 
convincing, in particular, since no parameter fitting was 
tried. 

Summarizing, we find convincing agreement with ex- 
periment for A as function of Ut as well as for / as a 
function of l/TZs- For the last, the available experimen- 
tal results allow to identify an almost quantitative agree- 
ment. The sensitivity of the experimental results on A as 
a function of l/TZs does not allow quantitative compari- 
son, and our results for / as a function of Ut deviate in 
their functional form from the available statements about 
experimental results. 



E. Mechanism of the oscillations, reaction-diffusion 
models and surface charge 

The voltage profiles U{t) in Figs. 2 and 5 show that 
there are two processes involved in the oscillations. 

The first process occurs on the slow time scale Ts of the 
semiconductor. It describes the exponential decay of the 
voltage Ut—U{T) — TZsj oc e~^^^" over the semiconductor 
layer according to Eq. (|23() . as long as the contribution 
of TZsj does not vary substantially. The decay time Tg is 
the Maxwell time due to resistance and capacitance of the 
semiconductor layer, accounts for the slow rise of the 
voltage U (t) over the gas discharge layer to a value above 
the current voltage characteristics of the gas discharge. 

The other process is the electric breakdown of the gas 
discharge layer for sufficiently large U{t) which leads to 
a current pulse and a rapid subsequent decay of U(^ 

It has been suggested by a number of authors 
17, 32, 33, 34, 35, 36, 32j that the current could be ap- 
proximated by a similarly simple equation of the type 
drj = g{U,j), where g vanishes on the current-voltage- 
characteristics. This would bring the equations into a 
reaction diffusion form. However, as we already have dis- 
cussed in |l2j . such an approximation of the underlying 
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equations lfTH|) - (j?T|) . and is not possible, since 
it would not admit the period doubling events observed 
in [l^ . and it would not allow the phase space plots in 
Figs. 3, 4 and 6 to intersect the characteristics with a 
nonvanishing derivative, as they definitely do. 

The physical reason for this behavior is the finite re- 
sponse time of the gas discharge layer, its "inertia" which 
doesn't allow an instantaneous reaction of the current. If 
ions are created by bulk impact ionization close to the 
anode, they will cross the whole gap until they reach 
the cathode and possibly liberate more electrons by sec- 
ondary emission. The time that the ions need to cross 
the gap, is therefore an important scale of internal mem- 
ory of the gas discharge. It can be approximated as 
Tion ~ L/{p.£) w L'^ / {plU{t)) where \£\ is some aver- 
age field within the gas gap. For the gap of L = 36 
[d = \ mm), the ion crossing time is estimated as 2.6-10'' 
for U{t) = 14 or as 1.5 • 10"' for U{t) = 24 (which corre- 
sponds to 0.6 or 1 /.ts in dimensional units). This time is 
of the same order or larger than the duration of a current 
pulse, both in our numerical solutions and in the experi- 
mental results of Fig. 5 in [id. (For the experiments on 
the 0.5 mm gap of Fig. 4 in [l^, the situation seems to 
be different.) 

Finally, it has been suggested in [sl] that the surface 
charge on the interface between gas and semiconductor 
could play an important role, in a similar way as in AC 
discharges. This is certainly true, but the surface charge 
q{T) is not an independent variable. Rather it is fully 
determined by the solution discussed above through 

,(r)=e.^i-i^-£(i,r). (35) 

The assumption that this surface charge is the only rel- 
evant charge in the whole system doesn't lead to a sat- 
isfactory description either, but the space charges in the 
gas discharge layer have to be taken into account, too. 

IV. STABILITY ANALYSIS: METHOD 

The direct numerical solution of the dynamical prob- 
lem is a time consuming procedure, that does not allow 
the exploration of a wide set of parameter values. We 
therefore have developed a linear stability analysis of the 
stationary state. It determines whether the stationary 
state is dynamically unstable and how small perturba- 
tions of such a state grow. In the present section, we 
present the method, and in the following one the results. 

A. Problem setting and stationary solutions 

The dynamical equations from section II. C are sum- 
marized as 

dr(7 = dzje+ jeMS) , je^CrS, (36) 

drS = j{T) + - tx£d,£, (37) 



TsdrU{T) = Ut~U{T)~nsj{T), (38) 

= d;,(t){z,T) +£{z,t) (39) 
with the boundary conditions 

dr£iO,T) = j(T)-Je(0,r), (40) 

dr£{L,T) ^ j{T)-^ML,r), (41) 
7 

0(L,r) = ^U{t) , 0(0, r) =0. (42) 

The stationary solutions form the starting point of the 
perturbation analysis. They solve the equations 

dzjeO = -jeo a{£o), (43) 

fj.£o dz£o = Jo - (1 + M)jeo, (44) 

dz^o = -£o, (45) 

Uo = Ut- n,jo (46) 

with boundary conditions 

Jeo(O) = JO , i^Jeo(i)=JO, (47) 

7 

MO) = , ML) = -Uo. (48) 



Eqs. lO-lESi with (gTI) and igHl define the current volt- 
age characteristics U = of a stationary discharge 
in the regime between Townsend and glow discharge 
.27, 28, , 30]. Eq. (gSJ is the load line due to the external 
circuit. The intersection of load line and characteristics 
defines a generically discrete number of stationary solu- 
tions of the system as a whole. 



B. Linear perturbations 

For linear perturbations about this stationary state, 
we use the ansatz 



jeiz,T) = Jeoiz) + Jel{z) (49) 

£iz,T) = £o{z) + £iiz) (50) 

(/.(z,t) = Mz) + Mz)e^\ (51) 

j{r) = Jo+Ji e^". (52) 



The lower index denotes the unperturbed stationary so- 
lutions while the lower index 1 denotes the linear pertur- 
bations about this stationary solution. The factorization 
of the perturbation into a z dependent function times the 
exponential e^''' anticipates the eigenvalue problem of the 
solution. 

In terms of the original variables, the explicit expan- 
sion in first order perturbation theory is a lengthy ex- 
pression, but in terms of the variables 

ai£o + ao£i jeijz) e^'' 

h = = — . - , and g = £o £i (53) 

cro£o JeO(Z) 
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the equations have a more compact form 



A 



dzQ = 
d.ji = 0, 



-h 



■^0 



.9 + - 



Ji 



with boundary conditions 

MO) = 0, 



A 



fo(0) 
A 



9(0)+ Jo h{0), 
g(L)+ JO h(L), 



(54) 

(55) 
(56) 
(57) 

(58) 
(59) 



(60) 
(61) 



(1 + At,) ML). 

Here the equation dzji — for the conservation of the 
total current is written exphcitly in order to bring the 
equations into the homogeneous form 



h 
9 



]eO 



V 



A 


_ J_ 



o\ 

1 


0/ 



9 
ji 

V01 



(62) 

The boundary conditions H58|l and H59|l at z = can 
be written as orthogonahty relations 



f jo \ / h 



£o(0) 

V / 



9 

h 







h 






9 


(?) 




ji 




I 


^1 



0. (63) 



The general solution v{z) of H62I) is therefore a superposi- 
tion of two independent solutions vi{z) and ^2(2) of H62|l 
that both obey in z = 0: 



I hiz) 
9{z) 
Ji{^) 

\Mz) 



= Ci Z7l(z) + C2 V2iz). 



(64) 



As initial conditions, one can choose, e.g., 


1 

\ 



MO) 



, V2{0) 



( \ 

£0(0) 

A 
1 



(65) 



V / 

The components of the two solutions are denoted as 
Vi{z) = (ft.i(0),5i(z),ji,i(z),0i^i(2;)). 

The boundary conditions (|60|l and (|61|l at z = L also 
have the form of orthogonality relations 




, 
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I ^'^ 


+ At, / 


\0i 



(66) 



Now each one of these two conditions determines the ratio 
C1/C2 of the general solution 

Ci \johi{L) + j^^g^{L) - jML) 

+C2 [3oh2{L) + j^g2{L)~jML)]^0, (67) 

Ci [-7^,il4(L) + (l + AT,)0l,l(L)] 

+C2 [-7^,Jl,2(i) + (l + AT,)0l,2(i)] = 0, (68) 

where ji^i(L) = 1 = Ji,2(i), since these components have 
this value at z = according to (jHSJ: ji,i(0) = 1 = 
Ji,2(0), and since the equation of motion for ji is dzji — 
0. A nontrivial solution of both (|67|l and (|68|) requires 
the determinant 

A = (69) 
:ioML) + £^)9i{L) - 1 ML) + £;;p:).92(i) - 1 

-Us + (1 + AT3)0i^i(i) -Us + (1 + AT,)0i,2(i) 

to vanish. This condition leads to a quadratic equation 
for the eigenvalue A. 



C. Rescaling with /i and numerical calculation 

The eigenvalue A can now be calculated numerically. 

First, it should be noted, that the equation of motion 
(I62|l has matrix elements of very different size, since ^ is 
a very small parameter. However, this apparent stiffness 
of the problem can be removed by introducing the new 
parameters 



_ _ J_ 'U 

A 

Ts =Ts^J, , S = -. 



(70) 



The introduction of rescaled current density and time 
scale and resistivity has a direct physical motivation. 
Previous analysis of the stationary solutions [23, 113, 113 
as well as the dynamical solutions of Section HI and [ij 
show that velocities should actually be measured on the 
time scale of the ions and not of the electrons. So the time 
scale should be measured in units of i+ = l/(ao^+i?o) — 
to/ ^ rather than in units of to = l/(aoMe^o)- The rescal- 
ing (|70|l directly follows from this consideration. 

Now the eigenvalue s can be calculated numerically as 
follows: First an initial estimate sq is chosen. Then the 
two initial conditions (|65|l at z = are integrated nu- 
merically with (|62|l up to z = L. Generically, the deter- 
minant A l|69|l will then be non-vanishing. The request 
that the determinant does vanish, fixes a new value for 
s that is used for the next step of the iteration within 
an under-relaxation method that garantuees the stabil- 
ity of the convergence. This procedure is repeated until 
an accuracy of 



Sk+l — Sk 
Sk+l 



(71) 
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is reached. 

The eigenvalue s is in general a complex parameter 
whose real part describes the growth or decay of the os- 
cillation amplitude while its imaginary part describes the 
oscillation frequency. Since s is a parameter in the equa- 
tion of motion (|62|l . also the vector v(z) has complex 
entries. Therefore 16 real functions Re hi{z), Im hi(z) 
etc. have to be integrated over z. It is convenient to 
also integrate the two real functions jeo and So that en- 
ter the matrix (|ti2|l together with the perturbations. The 
iteration program is written in fortran 90 with complex 
variables. For the integration of equations, a 4th order 
Runge-Kutta method is used. The number of grid points 
used was 500, since 1000 or 2000 grid points give essen- 
tially the same result. 



V. STABILITY ANALYSIS: RESULTS 

In the present section, the validity of the stability anal- 
ysis results are confirmed by comparison with numerical 
solutions of the full dynamical problem. The stability 
analysis is then used to determine the phase diagram for 
the onset of oscillating solutions. These phase diagrams 
are then compared with experimental results, again with 
semi-quantitative agreement. 



A. The structure of the results 

The stability analysis determines not only the complex 
eigenvalue A, but also the whole linear correction 



viz) = C\ 



vi{z) 



Us + {I + \t,) UigjL) 
TZs + [1 + Xn) Ui,i{L) 



V2{z) 



(72) 

up to the arbitrary complex constant Ci . 

This v{z) determines the evolution of current and volt- 
age in linear approximation about the stationary solution 
(jo,Wo): 



jir) = jo+ji e^"" + c.c. 
U(t) = Uo+Ui e^^ +C.C 



(73) 
(74) 



where c.c. denotes the complex conjugate. The ratio be- 
tween Ui and ji is fixed through the boundary condition 
(|()T|l to the value 



where r = 



Ji 



(1 + AT,)/7^, 

|l + Ar,| 
1 + 



(75) 



and 
Re At, 



(76) 



|1 + At, 



sma 



Im Ats 

11 + At, 



The final result is 



J(r) 
U{t) 



ere 



Rc At 
Ro At 



jo + e fi e^" cos(Im At 



cos(Im At 



- a + ao) 



(77) 
(78) 



where amplitude c and absolute phase ao reflect the arbi- 
trary factor Ci in (|72|l and are adjustable while all other 
parameters are fixed. 



B. Comparison with solutions of the full PDE's 

As a check of accuracy, these solutions are now first 
compared with numerical solutions of the full PDE prob- 
lem. 

For the set of parameters from Figs. 2, 3 and 4, the 
stationary solution is {joMo) — (1-49 • 10^^, 13.583), and 
the eigenvalue A has the complex value A = — 2.913 • 
10-6 ± i 4.822 • 10"^ As t^ = 340//i and 7^, = lAm/^i, 
the ratio of current and voltage amplitude r = 295//Lt 
and the phase shift a — 98.69° are determined through 
Eq. |7B1). 

The comparison of these predictions from the stability 
analysis with numerical solutions of the full PDE problem 
are shown in Fig. 8. Here the free parameters for the total 
amplitude c and the absolute phase ao were chosen such 
as to fit the PDE-data well. 



,x 10 




X 10 



FIG. 8: Comparison of j(T) and U{t)) from the stability anal- 
ysis (solid lines) with the result from the simulation (dashed 
lines) for the parameter values of Figs. 2-4. 



This visual agreement can be tested in more detail. In 
particular, we used the PDE-data in the time interval 
5 • 10^ < T < 6.5 • 10'"' to determine the phase shift a 
between Ui and ji. It is a = (100 ± 0.4)°, convincingly 
close to the predicted value of a = 98.69°. 

Increasing the total applied voltage Ut, the real part 
of the eigenvalue A grows until it becomes positive. 
This means that the stationary solution becomes lin- 
early unstable and perturbations will grow. An ex- 
ample of such behavior occurs for Ut = 24 with all 
other parameters as before. The stationary solution is 
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= (2.64 ■ 10-5,13.441), the eigenvalue is 
^ ±i 7.375 • IQ-^, the ratio of current and 



192//i and the phase shift is 



then (jo, Wo) 
A = 2.493-10 
voltage amplitude is r 
a = 99.83°. 

Fig. 9 shows again the comparison between these re- 
sults and the numerical solutions of the full PDE's 
Again, the agreement is very convincing. 
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FIG. 9: Comparison of j(r) and W(r)) from the stability anab 
ysis (solid lines) with the result from the simulation (dashed 
lines) for the parameter values of Figs. 5-6 where the station- 
ary solution is unstable. 



Of course, the predictive power of linear stability anal- 
ysis is limited to small perturbations with ji <ti jo and 
Ui <^ Uq. When the amplitude of the oscillation from 
Fig. 6 increases further, nonlinear couplings set in and 
the system finally reaches a limit cycle as shown in Fig. 7. 



C. Calculation of phase diagrams 

The stability analysis now allows one to derive the 
bifurcation line where a homogeneous stationary state 
looses its stability. Fig. 10 shows this bifurcation line 
for the parameters as a function of applied voltage 
Ut and conductivity 1/TZs for three different values of 7. 
Besides the value 7 — 0.08 used everywhere else in the 
paper, also results for 7 = 0.04 and 0.16 are shown to 
illustrate the sensitivity of theoretical predictions to this 
parameter. For Re A < 0, the stationary state is linearly 
stable, while for Re A > 0, the system is always in the 
oscillating state. 

Comparison with the experimental phase diagram in 
Fig. 11 for the gas gap with a corresponding width of 
d = 1 mm [TH l39| shows qualitative and quantitative 
correspondences, but also deviations. Experiments in the 
1 mm gap for Ut < 585 V {Ut = 20.5) do not exhibit os- 
cillations. The same holds theoretically for a secondary 
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FIG. 10: Bifurcation diagram for the parameters from I28II 
(where L = 36) and 3 different values of 7. The lines separate 
regions with Re A < where the stationary state is linearly 
stable from regions with Re A > where the homogeneous 
stationary state looses its stability. 
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FIG. 11: Blow up of the bifurcation diagram in Fig. 1101 for 
two different values of 7 and comparison with experimental 
data nil. I39t . Theoretical lines and experimental lines are in 
same region of parameters and have same limits. 



emission coefficient of 0.08 or smaller. In detail, experi- 
ments show that the raising phase transition line initially 
raises with positive slope then changes gradually to be- 
ing almost parallel to the as axis and then continues with 
negative slope up to the maximal experimentally reached 

For the low conductivity of the semiconductor layer, 
the experiment shows another bifurcation line almost 
parallel to the Ut axis at values of cr^ around 0.5 ■ 
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10~^/(ri cm). In dimensionless units this corresponds 
to a plateau at values of 1/7^^ around 0.5 • 10~^. An 
approach to such a plateau can also be seen in the cal- 
culated phase diagram. However, the theoretical curve 
crosses over continuously to this plateau, while the ex- 
perimental curve seems to show the intersection of two 
bifurcation lines with quite distinct slope. We have no 
explanation for this deviation. 

It is remarkable that the bifurcation theory also covers 
the almost horizontal bifurcation line for small l/TZg. An- 
other explanation for this experimentally observed fea- 
ture of the phase diagram would have been a breakdown 
of the continuum approximation: the recovery phase of 
the oscillation would have carried such a low current that 
the discreteness of the electrons would have to be taken 
into account. 

Finally, it was observed experimentally [Tol fllj that 
increasing the system size L while keeping other condi- 
tions unchanged, the frequency decreases and oscillations 
set in at higher voltages. This agrees with our calculated 
phase diagram in Fig. 11. Indeed, for Ut < 22.5, the 
homogeneous stationary state is stable for L = 72. 
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FIG. 12: Bifurcation diagram for 7 — 0.08, L = 36 and L 
72. 



VI. CONCLUSION 



We have analyzed the simplest model for a one- 
dimensional short gas discharge coupled to an external 
circuit with resistor, capacitance and stationary voltage. 
This analysis is directly applicable to experiments per- 
formed in jTol[lll| . 

We have presented fully numerical solutions as well as 
a linear stability analysis of the stationary state of the 
system which are in very good mutual agreement. The 
numerical solutions reproduce experimental observations 
of bistability and oscillations in a semi-quantitative man- 
ner, though the model is minimal and no attempt of pa- 
rameter fitting has been made. The stability analysis 
allows us to derive bifurcation diagrams in a simple man- 
ner, they also agree overall with experimentally obtained 
bifurcation diagrams. 

It should be remarked that we have constrained the 
analysis to the gap of 1 mm wide; the gap of 0.5 mm 
is so sensitive to the actual value of secondary emission 
7 that quantitative analysis based on a fixed value of 7 
seemed doubtful. 

We have reproduced a number of experimental obser- 
vations up to the dependence of oscillation amplitude on 
applied potential and of the oscillation frequency on the 
conductivity of the semiconductor layer, while discrep- 
ancies of other observables will stay a subject of inves- 
tigation. This opens up the way to investigate now the 
spatial and spatio-temporal patterns in the next step. 
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